
  debug = True 

  filenama = "~/input.nc" 
  filenamb = "~/st.nc" 
  filenamo = "out.nc" 

  system("after "+filenama+" "+filenamb+" < nml_T") 

  fila = addfile(filenama,"r") 
  filb = addfile(filenamb,"r") 

  hyam = fila->hyam 
  hybm = fila->hybm 
  ps   = fila->aps
  v4d  = fila->aclcac
  st   = filb->st

  nlev = dimsizes(hyam) 
  ntim = dimsizes(ps(:,0,0)) 
  nlat = dimsizes(ps(0,:,0)) 
  nlon = dimsizes(ps(0,0,:)) 

  pres = v4d
  pres = 0.
 
  rho  = v4d
  rho  = 0. 

  print("")  
  print("")  
  print("")  
  do k = 0, nlev-1
    pres(:,k,:,:) = ps(:,:,:) * doubletofloat(hybm(k)) + doubletofloat(hyam(k)) 
    if (debug) then 
        print("MAX,MIN,AVG of Pressure = "+max(pres(:,k,:,:))+"   "+min(pres(:,k,:,:))+"   "+avg(pres(:,k,:,:)))
    end if
  end do
  print("")  
  print("")  
  print("")  

  
  R = 8.3143

  rho = pres * 29. / ( R * st * 1000. )    ; kg / m3

  system("rm "+filenamo) 

  filo = addfile(filenamo,"c") 

  filo->rho=rho 

